Improved determination of the WIMP mass from direct detection data 
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Direct detection experiments searching for weakly interacting massive particle (WIMP) dark mat- 
ter typically use a simplified model of the Galactic halo to derive parameter constraints. However, 
there is strong evidence that this Standard Halo Model is not a good approximation to our Galaxy. 
We discuss previous attempts to extract the WIMP mass, cross-section and speed distribution from 
direct detection data and show that these lead to significant biases in the reconstructed parameter 
values. We develop and test an alternative model-independent method based on parametrising the 
momentum distribution of the WIMPs. This allows us to limit the analysis only to those regions 
of momentum space to which the experiments are sensitive. The method can be applied to a single 
experiment to extract the maximum information from a dataset, encoding combined information on 
the degenerate WIMP mass and interaction cross-section in a single parameter. This degeneracy 
can be broken by including data from additional experiments, meaning that the WIMP mass and 
speed distribution can be recovered. We test the momentum parametrisation method using mock 
datasets from proposed ton-scale direct detection experiments, showing that it exhibits improved 
coverage properties over previous methods, as well as significantly reduced bias. We are also able 
to accurately reconstruct the shape of the WIMP speed distribution but distinguishing between 
different underlying distributions remains difficult. 

PACS numbers: 07.05.Kf,14.80.-j,95.35.+d,98.62.Gq 



I. INTRODUCTION 

A large abundance of particle dark matter (DM) is 
required in our Universe to explain many observations on 
scales from the galactic to the cosmological (for a review, 
see e.g. Ref. Thus far, these particles have only been 
observed via their gravitational interactions and there 
are many experiments - completed, in progress and in 
development - which aim to directly detect Galactic dark 
matter through its non-gravitational interactions in the 
laboratory. While there is no lack of well-motivated DM 
candidates (e.g. Ref. 0-0), the nature of the dark matter 
particles remains an open question. Different classes of 
experiment are sensitive to different candidates and it 
is necessary to rely on a wide range of experiments to 
constrain the exact properties of the dark matter. 

We focus here on the search for weakly interacting mas- 
sive particle (WIMP) dark matter. WIMP direct de- 
tection experiments aim to identify and measure the en- 
ergies of nuclear recoils induced by WIMP-nucleon col- 
lisions @, H| ■ Excesses above the expected backgrounds 
have been observed in the CoGeNT [|[lo[ and CRESST- 
II [Tl| experiments. In addition, an annual modulation 
in the event rate has been observed at 8.9a significance 
by the DAM A/LIBRA collaboration While these 

tentative signals have been cited as evidence for a light 
WIMP of mass ~ 10 GeV, the various results are diffi- 
cult to reconcile in simple models of dark matter (T3I — 
[l5j and are in significant tension with the null results 
presented by CDMS and XENON100 [13 ■ It is im- 
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portant to accurately extract the properties of the WIMP 
from such experiments in order to check consistency with 
other detection channels (such as indirect detection [l8| 
and constraints from colliders [IH), as well as to provide 
constraints on the underlying models of supcrsymmetry 

WIMP direct detection experiments are typically anal- 
ysed within the framework of the Standard Halo Model 
(SHM) in which the velocity of galactic dark matter par- 
ticles is assumed to follow a Maxwell-Boltzmann distri- 
bution. This is a highly simplified assumption and can 
lead to a bias in the values of the WIMP mass and 
cross-section measured in direct detection experiments 
(see Ref. [HHHi and references therein). 

There have been several attempts to develop model in- 
dependent analysis schemes in order to avoid this bias as 
well as directly measuring the WIMP velocity distribu- 
tion. Frandsen et al. [3, [22| attempt to factorise out the 
astrophysical uncertainties, allowing a upper bound to be 
placed on the WIMP cross-section, while Alves et al. [23| 
parametrise the velocity distribution in terms of Legen- 
dre polynomials and attempt to reconstruct the polyno- 
mial coefficients from directional experiments. However, 
these techniques must be applied at a given WIMP mass, 
which would need to be measured independently (for ex- 
ample, in a collider experiment). Drees & Shan |24l. [25| 
have proposed calculating moments of the velocity dis- 
tribution function directly from data and using this to 
obtain the WIMP mass. However, the inclusion of finite 
energy windows in the mock experiments leads to sig- 
nificant bias and unreliable error estimates. Strigari & 
Trotta [2(| suggest using a simplified model of the Milky 
Way halo and simultaneously fitting the WIMP and halo 
parameters. However, the effectiveness of this method 
depends on how closely the actual Galactic halo can be 
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described by such a model. 

Here, we dev elop further the methods of Peter [13, [28| 
(see also Ref. [231), who has attempted to reconstruct 
the WIMP mass, cross-section and distribution function 
simultaneously by using a simple empirical parametrisa- 
tion of the velocity distribution. The analysis of mock 
datasets using this method indicates a significant bias 
in the reconstructed parameters even for simple veloc- 
ity distributions. We present an alternative method 
based on parametrising the momentum distribution of 
the WIMPs, which allows us to probe only that region 
of parameter space to which the experiment is sensitive. 
We test the method using mock datasets from proposed 
ton-scale direct detection experiments to assess its effec- 
tiveness. 

Section [IT] of this paper presents the current framework 
under which direct detection data is analysed, highlight- 
ing the need for alternative methods. Sec. IIIII details 
the Markov Chain Monte Carlo techniques which will 
be used to assess the proposed model-independent meth- 
ods, as well as the benchmark parameters used to gen- 
erate mock datasets. The statistical properties of previ- 
ous speed paramctrisation methods will be explored in 
Sec. IIVI We then present the new proposed momentum 
paramctrisation method and how it can be applied to 
a single experiment (Sec. |VJ) and to combined datasets 
from several experiments (Sec. IVI[) . The results of large 
numbers of parameter reconstructions are also presented 
in these Sections, before moving on to summarise the 
significance of this study in Sec. IVII1 



II. DIRECT DETECTION ANALYSIS 

The differential event rate per unit detector mass for 
nuclear recoils of energy En due to WIMP elastic scat- 
tering is given by poj : 
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where po is the local DM mass density, m x is the 
WIMP mass, <j p is the DM-proton cross-section (at zero- 
momentum transfer) and \x xp is the WIMP-proton re- 
duced mass, /j, xp = (m x m p ) / (m x + m p ). Here, we have 
considered only the spin-independent contribution to the 
rate, which dominates for targets with large atomic mass 
number, A. Assuming that the coupling to protons and 
neutrons is identical, this leads to an A 2 enhancement in 
the scattering rate (20[ . Finally, the loss of coherence due 
to the finite size of the nucleus is accounted for by the 
form factor, F 2 (Er), which we take to have the Helm 
form [U- 

The local velocity distribution of galactic WIMPs en- 
ters into the event rate through r](v m i n ), the mean inverse 
speed: 



d J v, 
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where /(v) is the velocity distribution in the Earth 
frame and u m in is the minimum WIMP speed required 
to produce a recoil of energy Er. For a detector of 
nuclear mass, m^, and reduced WIMP-nucleon mass, 
^xN = ( m x ra iv) / ( TO x + 171 n)i this minimum required 
speed is 
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Each experiment is then sensitive to a particular range 
of WIMP speeds, v £ [v m - m (E m i n ), v min (E max )], where 
E min and E max are the lower and upper limits of the 
signal window respectively. WIMPs with speeds below 
Wmin(-E'min) do not have sufficient kinetic energy to induce 
nuclear recoils above the threshold energy of the detector. 

Data from direct detection experiments are typically 
analysed using the Standard Halo Model (SHM), in which 
the galactic halo is modelled as an isothermal sphere with 
Maxwellian velocity distribution: 
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A typical value for the velocity dispersion is a w 
156 km s _1 , in order to match the peak speed to the local 
circular speed: vq = y/2a ~ 220 km s _1 . The motion of 
the Earth with respect to the Galactic halo is encoded 
in vi ag , which incorporates the Sun's motion within the 
Galaxy and the Earth's orbital motion about the Sun. 
Here, we neglect the time variation of vi ag and use an 
average value of |vi ag | « 230km s _1 [20| . 

We consider only non-directional detectors, for which 
it is useful to define the dircctionally averaged local speed 
distribution, /(«): 



(5) 



We also define the '3-dimensional' (3-D) speed distribu- 
tion, f3(v) = v 2 f(v), such that fs(v) dv is the fraction of 
WIMPs with speeds between v and v + dv. 

Despite the Standard Halo Model's frequent use, there 
is strong evidence that the true distribution function of 
galactic WIMPs is non-Maxwcllian. Even assuming the 
Standard Halo Model, there remains a great deal of un- 
certainty in the correct parameter values which should 
be used (see e.g. [HI, H2]). In addition, high resolu- 
tion N-body simulations suggest that there can be sig- 
nificant deviations from the maxwellian velocity distri- 
bution [33|-[35|], including features in /(v) especially at 
high velocity. Some simulations also show evidence for a 
dark disk, corotating with the baryonic disk of the Milky 
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Way, caused by the tidal distruption of satellites which 
are preferentially dragged into the disk plane [H, [37j . 
It is therefore necessary to consider the presence of addi- 
tional velocity structure beyond the SHM and the impact 
this may have on parameter reconstruction. 

III. PARAMETER ESTIMATION METHOD 

We now present the benchmark experimental and the- 
oretical parameters which will be used to test various re- 
construction methods as alternatives to the simple SHM 
assumption. The Markov Chain Monte Carlo methods 
used to explore the posterior distribution of the parame- 
ters arc also introduced. 



A. Benchmark Parameters and Experiments 

As noted in Ref. (25[, it is impossible to estimate the 
WIMP mass from a single experiment if no assump- 
tions arc made about f(v), so we consider three next- 
generation detectors, modelled on experiments which are 
currently in development. Each experiment is charac- 
terised by a single (suitably averaged) target nucleus 
mass, mjv, a total detector mass, mdet, an effective ex- 
posure time, £ oxp , and a pair of energies, -E m i n and E max , 
which mark the extent of the signal region. We focus on 
a particular set of experimental parameters in order to 
provide a concrete example of how the WIMP parameters 
can be estimated accurately. Table Q] shows the experi- 
mental parameters used in this work, which are chosen 
to approximately match those used by Peter [28| . 

We assume perfect and uniform detection efficiency - 
that is, all signal events and no background events sur- 
vive analysis cuts. We also assume perfect energy res- 
olution. For a real experiment, these assumptions will 
almost certainly not hold, for example due to variations 
in the relative scintillation efficiency of Xenon [4l| , but 
the results presented here should be viewed as a proof of 
principle in the ideal case. 

Fig. Q] shows the minimum and maximum accessible 
WIMP speeds for each experiment. All three experi- 
ments rapidly become insensitive to WIMPs with speeds 
less than ~ 1000 km s _1 as the WIMP mass drops be- 
low m x ~ 10 GeV. According to the RAVE survey 
the local galactic escape speed lies in the interval 
v esc G [498,608] km s" 1 . The fastest WIMPs impinging 
on an experiment will then have speed v = v\ ag + v esc . 
This is consistent with results from N-body simulations 
which indicate local escape velocities (in the Earth frame) 
of around 800 km s _1 [33]. This suggests that the exper- 
iments considered here generically have a low sensitivity 
to such light WIMPs, producing too few events for accu- 
rate parameter reconstruction. 

We therefore use benchmark masses of 50 GeV and 
100 GeV. For WIMP masses heavier than this, u m i n be- 
comes roughly independent of m x , leading to large de- 
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FIG. 1. Range of accessible WIMP speeds for each of 
the three mock experiments: XENONlT-like (solid blue), 
SuperCDMS-like (dashed green) and WArP-like (dot-dashed 
red). Each pair of lines corresponds to the maximum and min- 
imum accessible WIMP speeds for a given experiment. The 
outermost dotted red lines show the accessible speeds for the 
adjusted parametrisation range described in Sec. IVII 

generacy in the m x — a p plane [111 0, |44|. This will 
mean that for very heavy WIMPs, it will only be pos- 
sible to place lower bounds on the WIMP mass. In 
fact, at the statistical precision presented here this is 
also true for most mock datasets generated with m x = 
100 GeV. As a benchmark cross-section, we use the value 
(T p = 10~ 44 cm 2 , which is still a factor of several be- 
low current exclusion limits. We assume that the local 
dark matter density is known exactly and use the value 
po = 0.3 GeV cm~ 3 . As will be explained in Sec. IVII the 
precise values of a p and po are not particularly impor- 
tant due to the degeneracy between these two parame- 
ters. The total number of events from all three detectors 
combined typically ranges from around 300 to 600 for the 
benchmark parameters considered here. 

In order to ensure the robustness of the method, we use 
three benchmark models for the velocity distribution: 

(i) the Standard Halo Model (SHM), with a = 
156 km s _1 and wi ag = 230 km s ; 

(ii) a 50% Standard Halo Model with a 50% contribu- 
tion from a dark disk (DD); 

(hi) Via Lactea II data (VL-2). 

We model the dark disk velocity distribution as a 
Maxwcllian with a = 50 km s _1 and ui ag = 60 km s _ , 
similar to the typical values obtained by Ref. (4f| . A 50% 
contribution from the dark disk is at the upper limit of 
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XEN0N1T 



SuperCDMS [39J 



WArP [40] 



Detector Target 
Nuclear Mass, m^/a,mu 
Detector Mass, mdot/kg 
Exposure Time, t exp /days 
Energy Range, [i? m in, E 

max ]/keV 



Xe 

131 

1000 

365 

[2,30] 



Ge 

73 

100 

109.5 

[10,100] 



Ar 
40 

1000 

365 

[30,130] 



TABLE I. Parameter values for the three mock experiments used in this work, chosen to closely match those used in Ref. 
The meanings of the experimental parameters are described in Sec. IIII Al 



the range presented by Ref. [37J and we consider this as 
an extreme case. The third benchmark is the distribution 
function as extracted from the Via Lactea 2 N-body sim- 
ulation [!(| and presented in Ref. [33| . It is averaged over 
galactic radius in the range 7.5 < R < 9.5 kpc and mea- 
sured in bins of width 10 m s _1 . These benchmark veloc- 
ity distributions are illustrated in Fig. [2j The Standard 
Halo Model has the flattest velocity distribution with a 
tail extending up to roughly 800 km s _1 . This leads to a 
flatter spectrum and a larger number of events at higher 
energies than for the other two benchmark models. The 
VL-2 distribution produces fewer events compared to the 
SHM, as fewer WIMPs have sufficient energy to overcome 
the energy thresholds of the detectors. Due to the much 
smaller value of v\ ag for the dark disk model, WIMPs 
typically have much lower speeds, leading to even fewer 
events and a steeper spectrum. 



B. Parameter Exploration 

We generate mock datasets for each set of benchmark 
WIMP parameters and then analyse these using Markov 
Chain Monte Carlo (MCMC) techniques. For a given 
set of physical parameters, #tme = { m xi u pi f( v )}- WG 
calculate the expected number of events for each detector, 
N c , to an accuracy of 0.1 events: 
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For each experiment, the number of events actually 
observed, N , is drawn from a Poisson distribution with 
mean N e . Energies for the N Q events are then drawn 
from the normalised differential event rate: 
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FIG. 2. 1-D and 3-D speed distributions, f(v) and f-i(v), and 
mean inverse speed, rj(v), for the 3 benchmark speed models 
with parameter values as given in Sec. IIII Al : Standard Halo 
Model (SHM - solid blue), Standard Halo Model + Dark Disk 
(DD - dashed green) and Via Lactea 2 (VL-2 - dotted red). 
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We then use the publicly available CosmoMC code 
[47j to generate Markov chains to explore the poste- 
rior distribution of the parameter space. This parameter 
space typically comprises m x , <r p and additional parame- 
ters associated with the speed or momentum distribution 



of the WIMPs. The likelihood function used to generate 
the Markov chains is the same form used by CDMS [48| 
and XENON100 [49[, which for a single experiment is: 
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The full likelihood, £, is then the product of the likeli- 
hoods for the three separate experiments. 

For each set of benchmark parameters, we generate 250 
mock datasets. For each realisation, we perform 3 x 10 5 
Markov chain steps using the Metropolis-Hasting algo- 
rithm. Each chain is then thinned by a factor of 50 
and the initial 10 4 chain positions removed as burn-in. 
The results obtained were robust to doubling of either 
the thinning factor or the number of positions taken as 
burn-in. For a well-mixed chain, the distribution of chain 
positions as a function of the parameters, n{&), should 
be proportional to the local likelihood: n(9) cx C{9). 
For flat priors, this is then proportional to the poste- 
rior probability of the parameters 9. However, in order 
to ensure that the posterior is well-explored (particularly 
in cases where it may be multimodal), we perform the 
MCMC at a temperature T = 2, using the high tempera- 
ture likelihood function Ct = C 1 ^ (501 . The distribution 
of chain positions obtained at the higher temperature is 
then ut (0) oc C X I T {9). We recover the distribution of 
positions at T = 1 by 'cooling' the chain: 



where 9 Tec and #truc are the reconstructed and true values 
of the parameter of interest respectively. The A statistic 
quantifies the statistical deviation of the reconstruction 
from the true value. For a large number of mock datasets, 
the distribution of A should have a mean of zero and a 
standard deviation of unity. 

We sample the WIMP mass and cross-section logarith- 
mically in the ranges [10, 1000] GeV and [10, 10000] x 
10~ 47 cm 2 respectively, with log-flat priors on both. Un- 
less otherwise stated, parameters associated with the 
WIMP distribution are sampled linearly between zero 
and their maximum possible value allowed by normali- 
sation, with linearly-flat priors. 



IV. SPEED PARAMETRISATION 

We follow the method of Peter [28| and parametrise 
the WIMP speed distribution (in the Earth frame), as a 
series of N bins of constant value, with bin edges {v-i}: 



n(9)^n T (9)C 1 - 1/T (9). 



N 



(9) 



In order to obtain parameter limits from the mock 
datasets, we construct minimal credible intervals for the 
parameters of interest [5l| . In order to construct a p% 
confidence interval, we find a value of the marginalised 
likelihood for which p% of chain samples exceed this 
value. All parameter values for which the marginalised 
likelihood exceeds this limit are included in the confi- 
dence interval. This allows us to investigate the coverage 
properties of different methods and parameters. For each 
mock dataset, we record whether or not the true parame- 
ter value lies within this confidence interval. If a method 
produces consistent error estimates, the true parameter 
value should lie within the p% confidence interval for p% 
of realisations. This is referred to as exact coverage. 

We note that MCMC is typically poor at finding the 
best-fit point in parameter space [47[, so in order to ob- 
tain an estimate of parameters, we use the mean likeli- 
hood, as described in Ref. [47|. We bin the parameter 
of interest and average the likelihood of all points within 
each bin. This mean likelihood is then smoothed on the 
width of the bins and the parameter value which max- 
imises the mean likelihood is taken as a best estimate of 
the underlying parameter. 

As an estimate of an error on a parameter, we calculate 
its standard deviation within a chain, a. While this does 
not correspond to the Gaussian error on the parameter 
(particularly in the case of multimodal and skewed pos- 
terior distributions), it can be used as a rough estimate 
of the uncertainty of a particular reconstruction. We can 
also use this to calculate the 'pull' statistic A (see e.g. 
Ref. H3), 
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where the top-hat function, W, is defined as: 
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1 v € [vi,Vi + Av] 
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Here, we consider N — 5 bins, and take the bins to be of 
constant width, Av = 200 km s _1 , covering the range v 
— 1000 km s _1 . Imposing normalisation of the speed 
distribution, we obtain the following constraint on the 
{9iY- 



N 



5> [(v i + Av) 3 -vf]/3 = l. 



(13) 



For notational convenience, we also define 

9i = 9l [(vi + Av) 3 - vf] /3 , 
such that the normalisation condition becomes 

N 



(14) 



(15) 



A. Results 

For comparison with later methods, we consider here 
a single benchmark model: m x = 50 GeV and the SHM. 
Fig. El shows the fitted values for the WIMP mass, m rcc , 
obtained from 250 mock datasets. This distribution 
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FIG. 3. WIMP masses reconstructed using the speed 
parametrisation method from 250 realisations. The bench- 
mark speed distribution is the SHM. The true mass of 50 
GeV is shown as a dashed vertical line. 



FIG. 4. Distribution of the A statistic, defined in Sec. IIIIBI 
for 250 realisations using the speed parametrisation method. 
The benchmark speed distribution is the SHM, with a 50 GeV 
WIMP. 



shows a peak around 45 GeV, as well as a significant num- 
ber of datasets reconstructed at ~ 100 GeV. As pointed 
out by Ref. HH , some mock datasets will not be represen- 
tative of the underlying benchmark parameters, having 
more events at high energies than expected, for exam- 
ple. This can lead to 'bad' reconstructions with a fitted 
WIMP mass higher than the benchmark value. Thus, the 
reconstructions near 100 GeV do not necessarily signify 
a failure of the reconstruction method. 

However, a coverage study of 68% and 95% confi- 
dence intervals for this method shows significant under- 
coverage: 36 ± 3% coverage and 63 ± 3% coverage respec- 
tively for the two intervals. This indicates that while the 
mass reconstructions appear to be distributed close to 
the true value, the corresponding error estimates must 
be too small. Fig. 2] shows the A distribution for the 250 
realisations of this method. The standard deviation of 
A is a - a = 3.5, which is significantly greater than unity, 
showing that the speed parametrisation method signifi- 
cantly underestimates the errors on reconstructed values. 
As demonstrated in Ref. [28| , the problem of poor recon- 
structions using this method does not appear to be signif- 
icantly improved by increasing the number of bins in the 
speed parametrisation and worsens for more complicated 
speed distributions. 

Fig. [5] shows the reconstructed speed distribution for a 
typical realisation using this method. The reconstructed 
mass value is log 10 (m rcc / GeV) = 1.48 ± 0.06, com- 
pared to the benchmark value log 10 (m x / GeV) = 1.699. 
The mean inverse speed is under-estimated in the range 
— 200 km s -1 and slightly over-estimated at higher 
speeds. However, the reduced m roc increases the mini- 
mum accessible speed of the experiments, meaning that 
the experiments are less sensitive to the shape of the 
speed distribution at low speeds. Moreover, a reduced 



value of the reconstructed mass serves to steepen the 
spectrum, reconciling the flattened ?7(vmin) at high speeds 
with the data. This is because varying the mass of 
the WIMP 'rescales' the spectrum, due to the relation 

E R OC MxA^min- 

In Fig. [6l we plot rj/m x as a function of recoil energy, 
En, for the SuperCDMS-likc experiment. We rescale rj 
by 1/ m x because this factor appears in the event rate 
and we are then able to compare the spectra of events 
from different models. The solid line shows the mean in- 
verse speed in the SHM, using the true WIMP mass of 50 
GeV. We also show a binned approximation to the SHM 
(dashed line) obtained using the true values of the bin 
parameters {gi} and the true WIMP mass. Finally, we 
show the reconstructed mean inverse speed (dot-dashed 
line) using the reconstructed WIMP mass of 30 GeV. We 
see that the binned approximation to the SHM, which 
should represent a 'good' reconstruction, actually recov- 
ers the spectrum poorly compared to the reconstructed 
values. In particular, we note the energy range of the ex- 
periment spans two bins in the binned approximation to 
the SHM, but three bins in the MCMC reconstruction, 
allowing a closer approximation to the true spectrum. 

Thus, the speed distribution parameters can be ex- 
plored to provide a good fit to the data, with the recon- 
structed mass varying to compensate. As can be seen 
from Fig. [5J for a fixed bin width in velocity space, the 
size of bins in energy space can be reduced by moving to 
lower masses. This should allow a closer fit to the data 
and may explain why there appears to be a bias towards 
lower mass values. 



7 



1 



0.8; T 




200 400 600 800 1000 



v / km s 

FIG. 5. Reconstructed speed distribution, f(v), and mean 
inverse speed, rj(v), using the speed parametrisation method. 
The benchmark model used was a 50 GeV WIMP with an 
SHM speed distribution. The upper pane shows the underly- 
ing SHM speed distribution (solid blue) and the fitted values 
of the speed bin parameters (red points). The errors on the 
bin values are within-chain standard deviations as described 
in Sec. IIIIBI The lower pane shows the mean inverse speed 
corresponding to these fitted values (dashed red line) and the 
true mean inverse speed (solid blue). 
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FIG. 6. The rescaled mean inverse speed, r//m x , measured 
in the SuperCDMS-like experiment as a function of recoil en- 
ergy, _Er. The same mock dataset was used as in Fig. [S] The 
underlying Standard Halo Model distribution (solid blue) uses 
the true WIMP mass of 50 GeV, as does the binned approx- 
imation to the SHM (dashed red). The reconstructed mean 
inverse speed (dot-dashed black) uses the reconstructed value 
of 30 GeV. 



the reduced momentum distribution, /(p): 
f°° /(p) 1 

ViPmin) = / d 3 p = ^(Pmin/Mviv)- ( 19 ) 

The event rate can be rewritten as: 



V. MOMENTUM PARAMETRISATION FOR A 
SINGLE EXPERIMENT 

When considering the speed distribution of the 
WIMPs, we see that each experiment has a different 
range of sensitivity and that varying the WIMP mass 
changes this range. However, we can instead consider a 
'reduced WIMP-nuclcon momentum,' 



Pn — Mx^ 1 



(16) 



defined separately for each target nucleus. We now note 
that the accessible range in p N for each experiment is 
independent of the WIMP mass: 



Pmin(ER) = M x AfW m i n (i?fl) 



m N E R 



(17) 



We therefore rewrite the differential event rate in terms 
of the new momentum variable: 

dE R 2[i l xp m x 
where fj is the mean inverse momentum associated with 
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where we have defined 



D(a p ,m x ,m N ) 



(Tp(l x N 



(20) 



(21) 



which encodes all information about the WIMP mass and 
cross-section and controls the overall scale of the event 
rate. 

We can again define a directionally averaged mo- 
mentum distribution, f(p) = j ~(p / 'h x n) I ' P? x n > an< ^ 
parametrise this in terms of 5 constant bins, with bin 
values {hi}. We parametrise f(j>) only over the range 
of sensitivity of the experiment: p £ [PajP&L where 
Pa, 5 = Pmin(-Emin,roax)- However, wc still wish to im- 
pose some normalisation constraint on the momentum 
distribution parameters. Each experiment now probes 
a well-defined (but unknown) fraction of WIMPs, a^, 
given by 



a N 



f(p) P 2 dp- 



(22) 
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The momentum parameters arc therefore normalised ac- 
cording to 



N 



5> 



OtN 



(23) 



i=i 



where hi is defined analogously to fji in Eq.rj5] We absorb 
the unknown a at into D, such that the momentum dis- 
tribution parameters, {hi /aw}, are normalised to unity 
and 



D(a p ,m x ,m N ) = a N ^ 



(24) 



Finally, it is necessary to introduce a parameter A 
which models the contribution from WIMPs above the 
maximum energy of the signal region: 



.4 



M d 3 p . 

P 



(25) 



Because the precise form of f(p) above the upper energy 
threshold is undetermined by the experiment, the contri- 
bution of A to the normalisation, a at, cannot be calcu- 
lated and is therefore not considered. Instead, we include 
conservative constraints on A such that its contribution 
alone cannot exceed the normalisation of f(p): 



A < (Pmin(-Emax)r 



We also note that 



___] d 3 p <^_, 

P Pb 



(26) 



(27) 



and thus impose the following additional constraint on 
the paramaters: 



— lv(Pa) - v(pb)] < — • 
ajv Pb 



(28) 



We therefore perform parameter reconstructions using 
the parameters D, {hi/aN} and A/chn ■ Because the frac- 
tion of high momentum WIMPs is expected to be rela- 
tively low, we sample the parameter A logarithmically, 
with a log-flat prior. 



A. Results 

We consider again a single set of benchmark parame- 
ters, namely a 50 GeV WIMP with an SHM speed dis- 
tribution. We apply the momentum parametrisation to 
mock datasets from the WArP-like Argon experiment. 
The reconstructed D values, D rcc , are shown in Fig. [7] in 
units of 10 7 cm 2 kg~ 2 . In all reconstructions, the poste- 
rior distribution is unimodal, having separate parameters 
to describe the scale (D) and shape ({hi}) of the event 




D /(10 7 cm 2 kg' 2 ) 



FIG. 7. Reconstructed values for the scale parameter, D ICC , 
for the Argon experiment using the momentum parametrisa- 
tion method from 250 realisations. The benchmark speed 
distribution is the SHM. The value of Dtruc = 75.6 x 
10 7 cm 2 kg -2 is shown as a dashed vertical line. 



rate. The number of reconstructions is peaked at the cor- 
rect value, however, the distribution does not appear to 
be symmetric. In fact, the average reconstructed value is 
log 10 (D) rcc = 1.865 ±0.004, compared to the input value 
of log 10 (£))t rU o = 1-878. This represents a slight bias (of 
less than 1%) towards smaller values of log 10 (£>). 

However, this is smaller than the typical statistical un- 
certainty in a single reconstruction, which is ~ 4%. In 
addition, this method results in overcoverage of the true 
parameter, with values of 76 ± 2% and 98 ± 1% respec- 
tively for the 68% and 95% confidence intervals. This 
method therefore allows us to place reliable conservative 
estimates on the parameter D. 

We show in Fig. [8] the reconstructed momentum 
distribution and mean inverse momentum for a typi- 
cal realisation, for which the reconstructed D value is 
1.81^0 05- The underlying momentum distribution has 
been rescaled by l/a\ r to allow a comparison to the 
reconstructed values. We see that the the momentum 
distribution is well reconstructed and the mean inverse 
momentum is accurately recovered at low and high mo- 
menta. In the middle of the momentum range, however, 
jy(Pmin) exceeds the true value. Because only a single 
experiment is being used, the measured spectrum is par- 
ticularly susceptible to Poisson fluctuations. The mock 
dataset used here has a slight excess of events around 
-Er ps 60 keV, corresponding to pAr ~ 30 MeV, which 
may explain the reconstructed excess. 

In addition, this may be a consequence of the partic- 
ular parametrisation. The constant-bin parametrisation 
of f(p) leads to a parametrised ?7(pmin) which is con- 
cave downwards in each bin, while the underlying func- 
tion is strictly convex downwards in this region. Thus, 
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2.0 r 




FIG. 8. Reconstructed momentum distribution for a single 
Argon experiment using a benchmark of a 50 GeV WIMP and 
the SHM. The upper pane shows the SHM momentum distri- 
bution (solid blue) and reconstructed bin values (red points). 
Because the posterior is unimodal, we also display vertical 
errorbars showing the extent of the 68% confidence region for 
each bin. Note that these errors are strongly correlated. The 
lower pane shows the corresponding reconstructed mean in- 
verse momentum (dashed red) and the mean inverse momen- 
tum in the SHM (solid blue) . The underlying distribution has 
been rescaled by 1/<xax for comparison to the reconstructed 
values. 

7?(Pmin) tends to be slightly overestimated, leading the 
scale parameter D to be reduced to compensate for this. 
With datasets containing more events, the number of bins 
could be increased, in order to reduce this bias on D and 
maintain it at below the level of the statistical uncer- 
tainty. 



VI. MOMENTUM PARAMETRISATION FOR 
SEVERAL EXPERIMENTS 

The reduced momentum method allows us to extract 
the maximum amount of information from a single ex- 
periment, making no assumptions about the underlying 
velocity (or momentum) distribution. However, informa- 
tion about the mass and cross-section are encoded in the 
parameter, D, and cannot be extracted using a single 
experiment alone. 

Because a different momentum variable pjy can be de- 



fined for each experiment, it is necessary to choose a sin- 
gle experiment and parametrise the momentum distri- 
bution defined with respect to that experiment. It may 
be necessary to adjust the lower and upper limits of the 
parametrisation (beyond the values of -Emm and i? mM 
used in the experiment) to accommodate as much of the 
data as possible from all experiments. In the single exper- 
iment considered in Section [Vj the WIMP-Ar momentum 
was parametrised in the range pAr G [23.6, 49.2] MeV, to 
match the sensitivity of the Argon experiment. How- 
ever, as can be seen in Fig. [TJ this sensitivity win- 
dow does not match that of the other experiments. If 
we extend this interval, and parametrise in the range 
PAi G [3.6, 53.0] MeV, we can enclose the sensitivity re- 
gions of all three experiments as closely as possible, as 
shown by the dotted curves in Fig. [TJ We again use 5 
bins in momentum space, with an additional parameter 
to control a constant offset. 

In theory, any of the three experiments could have been 
chosen to define the momentum variable. However, some 
choices of experiment are less practical. For example, in 
order to use the XENONlT-like experiment, it would be 
necessary to parametrise the momentum over the range 
pxc G [11,162] MeV. This is because at high WIMP 
masses the remaining two experiments have maximum 
accessible speeds of ~ 500 km s _1 . This corresponds to 
very high values of the WIMP-Xe reduced momentum 
because of Xenon's comparatively higher mass. Because 
this parametrisation would cover such a large momentum 
range and because structures in the distribution function 
tend to be concentrated at low speeds, a very large num- 
ber of bins would be required to accurately model such 
structures. 

By comparison, using the WArP-like Argon experi- 
ment allows us to parametrise only as much of the mo- 
mentum space as required to accommodate data from all 
three experiments. In the speed parametrisation method, 
the WIMP mass could be varied to adjust the range of 
speeds accessible to the experiments. This reduces the 
sensitivity of the likelihood to some of the speed bin pa- 
rameters, meaning that these parameters can be adjusted 
with little effect on the likelihood value. This results 
in a spurious freedom in the remaining bin parameters, 
which can be varied to achieve a good fit to the data. As 
demonstrated in Sec. IIV1 this results in a bias towards 
lower values of the reconstructed WIMP mass in order 
to reduce the size of the bins in energy space. Using 
the momentum parametrisation method, we reduce this 
effect by ensuring that the likelihood is sensitive to all 
momentum bin parameters as much as possible, as the 
accessible speed range of the analysis tracks more closely 
the accessible speed range of the experiments. In addi- 
tion, for a fixed bin width in momentum space, the bin 
width in energy space is much less sensitive to the WIMP 
mass. 

Unfortunately, this method does not allow the WIMP- 
nucleon cross-section to be extracted; because the con- 
tributing WIMP fraction, a, is unknown, we can only 
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FIG. 9. WIMP masses reconstructed using the speed and mo- 
mentum parametrisation methods from 250 realisation. The 
benchmark speed distribution is the SHM. The true mass of 
50 GeV is shown as a dashed vertical line. 



obtain a lower bound. This is a fundamental limitation of 
any method which makes no assumptions about the un- 
derlying speed distribution. Without knowing the frac- 
tion of WIMPs with speeds within the signal window of 
the experiment, we cannot determine the cross-section. 
However, the cross-section appears in the event rate only 
through the degenerate combination a p po. Estimates of 
po from mass modelling of the Milky Way and from stel- 
lar kinematics tend to have a large uncertainty, with po 
topically lying in the range 0.2 — 0.4 GeV cm~ 3 (e.g. Ref. 
[HjjlHH). Thus, any estimate of the WIMP-nucleon inter- 
action cross-section would have an inherent uncertainty 
in any case. 



A. Results 

We first compare results for the momentum method 
to the speed parametrisation method presented in Sec- 
tion IIVI We use the same mock datasets generated for 
the 50 GeV, SHM benchmark presented previously. The 
results of both the momentum and speed methods are 
shown in Fig. |9] In the case of the momentum method, 
the distribution of realisations is now more closely peaked 
around the true mass of 50 GeV. Furthermore, the mo- 
mentum method produces substantially improved cover- 
age properties, as summarised in Table ITT1 

Fig. [TU] shows the reconstructed WIMP- Argon mo- 
mentum distribution using the same mock dataset as 
used for Fig. [5] The benchmark distributions have 
been rescaled by a so that they can be compared 
to the reconstructed values. In this case, a = 





Speed Method 


Momentum Method 


68% Coverage 


36 ± 3% 


71 ± 3% 


95% Coverage 


63 ± 3% 


92 ± 2% 



TABLE II. Coverage statistics for the speed and momen- 
tum parametrisation methods for a 50 GeV SHM benchmark 
model. 
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FIG. 10. Reconstructed momentum distribution from all 
three mock experiments using a benchmark of a 50 GeV 
WIMP and the SHM. The upper pane shows the SHM mo- 
mentum distribution (solid blue) and reconstructed bin val- 
ues (red points). The errors on the bin values are within- 
chain standard deviations as described in Sec. MI Bl The lower 
pane shows the corresponding reconstructed mean inverse mo- 
mentum (dashed red) and the mean inverse momentum in 
the SHM (solid blue). The reconstructed values have been 
rescaled by a for comparison to the true distribution. 



0.995, so we can reconstruct both the mass and cross 
section accurately: log 10 (m roc / GeV) = 1.62 ± 0.31 
and log 10 ((Tp/10~ 47 cm 2 ) = 2.99 ± 0.18, compared to 
the true values of log 10 (m tr ue/ GeV) = 1.699 and 
log 10 ((7p/10 -47 cm 2 ) = 3.0. While there is no way to 
know a priori whether a will be close to unity, the ac- 
curate reconstruction of the mass, cross-section and mo- 
mentum distribution show that momentum parametrisa- 
tion can offer a significant improvement over the speed 
parametrisation method. 

We now present the results of the momentum method 
for all six sets of benchmark parameters, 50 GeV and 100 
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FIG. 11. Distribution of reconstructed masses, m rcc , using the 
momentum method for 250 reconstructions. The true mass 
of 50 GeV is shown as a dashed vertical line. 



WIMP Mass 
50 GeV 100 GeV 



SHM 



DD 



VL-2 



71 ± 3% 
92 ± 2% 

61 ± 3% 
94 ± 1% 

61 ± 3% 
86 ± 2% 



65 ± 3% 
94 ± 1% 

58 ± 3% 

91 ± 2% 

62 ± 3% 

92 ± 2% 



TABLE III. 68% and 95% confidence interval coverage results 
for the momentum parametrisation method using a variety of 
benchmark parameters, as defined in Sec. MI Al 



GeV WIMPs, with Standard Halo Model, dark disk and 
Via Lactea-2 distribution functions. The distributions of 
reconstructed masses are shown in Fig.QTJfor the 50 GeV 
WIMP and Fig. [H for the 100 GeV WIMP. 

For the 50 GeV benchmark, the distribution of recon- 
structions is peaked at the true value, though in all three 
cases there are a number of datascts reconstructed at 



FIG. 12. As Fig. [TT] for m x = 100 GeV. 



higher masses. For some of the mock datasets, the pos- 
terior distribution for the WIMP mass is multimodal, 
with a peak near the true value as well as a peak above 
~ 100 GeV. For reconstructions using a fixed speed 
(or momentum) distribution, these would correspond to 
'bad' reconstructions, as mentioned previously, in which 
the spectrum of events is flatter than expected. When 
the momentum distribution is allowed to vary, as here, 
the event rate can be well fit by more than one region of 
the mass parameter space. 

For the 100 GeV benchmark, the SHM and VL-2 mod- 
els show similar structures, with a broad peak of re- 
constructions at or near the correct values, as well as 
a smaller tail up to masses of 1000 GeV, the upper limit 
of the prior. The 100 GeV datasets contain fewer events 
than their 50 GeV counterparts, so we would expect the 
spread of reconstructed values to be broader. Also, as 
previously noted, as the WIMP mass exceeds the mass 
of the target nucleus, the shape of the event spectrum 
becomes roughly independent of the WIMP mass. The 
largest nuclear mass used here is ^4xo = 131, mean- 
ing that for WIMP masses significantly above m x w 
131 amu »s 122 GeV, the underlying value is difficult to 
reconstruct exactly. For these benchmark parameters, 
we can therefore only place a lower bound on the WIMP 
mass and when calculating coverage statistics, we use 1- 
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tailed limits (i.e. ap% confidence limit encloses 
of the marginalised posterior). 

We report coverage statistics for the various bench- 
mark parameters in Table Mil For the Standard Halo 
Model there is approximately exact coverage for both 50 
and 100 GeV WIMPs. The remaining benchmark pa- 
rameters display some undercoverage, though still much 
improved over that achieved by the speed parametrisa- 
tion method. The poorest coverage is achieved for the 100 
GeV DD benchmark, for which the 68% confidence inter- 
val has a coverage of 58±3%. This is to be expected from 
the poorly distributed reconstructions shown in Fig. [T^l 
For the 100 GeV dark disk benchmark, there appears to 
be a significant bias in the distribution of reconstructed 
values, which peaks around 70 GeV. We explore the ori- 
gin of this bias in the next section, where we examine the 
speed distributions reconstructed using this method. 



B. Recovering the speed distribution 

We will now consider how the speed distribution can 
be reconstructed from the momentum paramctrisation. 
For a set of constant bins in momentum space, the po- 
sitions and widths of bins in velocity space is dependent 
on the WIMP mass. It is therefore difficult to extract 
precise statistical information on the speed distribution, 
as the bin values will be very strongly correlated with the 
WIMP mass. Instead, we take the reconstructed WIMP 
mass as fixed and use this to obtain a speed distribution 
from the momentum distribution parameters. Without 
treating the covariance of the WIMP mass and the bin 
parameters in full, the reconstructed speed distribution 
will depend strongly on the reconstructed mass value. 
However, this naive approach should give an indication 
of whether accurate reconstructions are possible. 

First, we consider a 50 GeV WIMP with SHM dis- 
tribution, as an archetypal WIMP model with a well- 
behaved distribution function. We show a typical recon- 
structed speed distribution in Fig. II 31 using the same 
mock dataset as Fig. [TUJ In this case, the reconstructed 
value of TO rec is 42 GeV and the speed distribution ap- 
pears to be accurately reconstructed within the error es- 
timates. 

Next, we consider a reconstruction for a 100 GeV 
WIMP with DD distribution function. One example is 
shown in the left-hand panes of Fig. [T4J for a dataset 
with reconstructed mass log 10 (TO rec /GeV) = 1.83 ± 0.15, 
compared to the true value of log 10 (m x /GeV) = 2. The 
speed distribution appears to be well recovered at all 
speeds. However, there is a significant discrepancy in the 
mean inverse speed below ~ 150 km s~ . This is because 
the DD distribution function is very rapidly varying at 
low v, meaning that the ansatz of constant bins can no 
longer be applied. As observed in the speed parametri- 
sation method, the event spectrum can be steepened by 
moving to lower mass values and this may explain why 
there is significant bias and poor coverage for this set of 
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FIG. 13. Reconstructed speed distribution from all three 
mock experiments using the momentum parametrisation 
method. The benchmark is a 50 GeV WIMP and the SHM 
distribution function. The upper pane shows the underlying 
SHM speed distribution (solid blue) and the fitted values of 
the speed bin parameters (red points). The errors on the bin 
values are within-chain standard deviations as described in 
Sec. IIII Bl The lower pane shows the mean inverse speed cor- 
responding to these fitted values (dashed red line) and the true 
mean inverse speed (solid blue). The underlying distributions 
have been rescaled by a for comparison to the reconstructions. 
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FIG. 14. As Fig. [TJ for a 100 GeV WIMP with DD distri- 
bution function using 5 momentum bins (left panes) and 7 
momentum bins (right panes). 
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FIG. 15. Distribution of reconstructed masses using the 7-bin 
momentum method for 250 reconstructions for a DD bench- 
mark distribution. The true mass of 100 GeV is shown as a 
dashed vertical line. 



benchmark parameters. 

In the right-hand panes of Fig. [HI we show results 
from the same mock dataset reconstructed using 7 bins 
in momentum space. The reconstructed mass is now 
log 10 (m roc /GeV) = 2.21 ± 0.27, with the mean inverse 
momentum more closely reconstructed than for the 5 bin 
case. Figure [T5] shows the distribution of reconstructed 
masses for a 100 GeV WIMP with a DD distribution 
function using 7 bins in momentum space. The recon- 
structed masses are now more broadly distributed around 
the benchmark value, with improved coverage compared 
to the 5 bin case: 67 ± 3% and 94 ± 1%. We have found 
that increasing the number of bins for the 50 GeV SHM 
benchmark leaves the coverage properties and distribu- 
tion of reconstructions largely unchanged, indicating that 
increasing the number of bins can be used to check the 
robustness of the reconstructions. 

Finally, we consider the discriminatory power of the 
reconstructions. Returning to the 50 GeV SHM bench- 
mark, we plot a single speed distribution reconstruction 
in Fig. 1161 as well as all three benchmark speed distribu- 
tions for comparison. The reconstruction is reasonably 
consistent with both the SHM and VL-2 models and dis- 
plays only mild tension with the DD model. In addition, 
the benchmark distributions in Fig. ll6l have been rescalcd 
by the true value of a for comparison with the recon- 
structed values. In a real experiment, the value of a is 
unknown, further reducing the potential to discriminate 
between different models. Only in the case of more ex- 
treme distribution functions, such as a dark disk, might 
it be possible to make a distinction between the many 
possible underlying models. Thus, while the momentum 
paramctrisation method can provide good constraints on 
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FIG. 16. Reconstructed speed distribution from all three 
mock experiments using a benchmark of a 50 GeV WIMP 
with SHM distribution. The reconstructed values have been 
rescaled by a for comparison to the true distribution. The 
three different benchmark speed distributions defined in Sec. 
IIII Al have been overlaid: SHM (solid blue), DD (dashed 
green) and VL-2 (dotted red). 



the mass of the WIMP, it remains difficult to probe the 
speed distribution function. 



VII. CONCLUSION 

We have presented a novel method of analysing direct 
detection datascts which aims to reconstruct as much 
information as possible about the WIMP mass, cross- 
section and distribution function while making no as- 
sumptions about the shape of the underlying speed dis- 
tribution. To do this, we parametrise the WIMP mo- 
mentum distribution using a simple empirical paramctri- 
sation. This ensures that we do not parametrise those 
regions of velocity space to which the experiments are 
not sensitive, thus preventing spurious reconstructions, 
as seen previously in the use of speed parametrisation 
methods. 

In the case of a single experiment, this method can be 
applied exactly and allows one to extract the maximum 
possible information, at the cost of losing access to in- 
formation about either the WIMP mass or cross-section 
separately The overall 'scale' parameter £>, defined in 
Eq. [M] which encodes information about both the mass 
and cross-section, can be reconstructed with only a small 
bias, and conservative limits can reliably be placed on D. 
The D values from many different experiments can then 
potentially be used to place bounds on the values of the 
WIMP mass and cross-section. The momentum distribu- 
tion is also accurately reconstructed, though an increased 
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number of parameter bins may be required as the number 
of signal events increases. 

For multiple experiments, the range of the paramctri- 
sation must be extended to cover the sensitivity regions of 
all experiments. For estimation of the WIMP mass, this 
allows us to achieve significant improvements in cover- 
age and reduction in bias over previous methods, though 
the momentum paramctrisation method still suffers from 
some under-coverage for more extreme distribution func- 
tions, such as the dark disk. In these cases, the coverage 
properties can be improved by using more bins in mo- 
mentum space. Without making any assumptions about 
the WIMP speed distribution, however, we cannot es- 
timate the interaction cross-section due to its degener- 
acy with the fraction of WIMPs accessible to the exper- 
iments. This is an unavoidable problem, but is rendered 
somewhat irrelevant by large uncertainties in the local 
DM density. 

Reconstruction of the WIMP speed distribution re- 



mains difficult. The finite sensitivity window of direct de- 
tection experiments means that information on the nor- 
malisation of f(v) is lost, making comparison to theo- 
retical models difficult. At the event rates studied here, 
it does not appear to be possible to distinguish between 
different distribution functions. If a clear dark matter 
signal is observed at next-generation ton-scale detectors, 
however, it should still be possible to accurately estimate 
the WIMP mass in a model-independent fashion. 
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